2. Species-area relationships in the Aegean Islands

Hammoud et al (2021) examined the relationship between island area and species richness in the Aegean Islands across multiple plant and animal taxa, taking into account historical sea level change over the Quaternary. They find that land-bridge islands, which were historically connected to the mainland, exhibit weaker correlations between island area and endemic species richness relative to true islands, which were never connected with the mainland. The results of this study suggest that periodic land bridge formation between islands and the mainland may inhibit the differentiation of insular populations. For the sake of this vignette, however, we will only be looking at species-area relationships for angiosperms on true islands in the Aegean archipelago.

First, we need to download bathymetry data for the Aegean Islands. For this vignette, we will use data from the Generalised Bathymetric Chart of the Oceans (GEBCO)

Now we can generate the interval file. For the sake of this demonstration, we will generate an interval file with sea level depth intervals (as opposed to time intervals). with a time cutoff of 20 kya (the approximate last glacial maximum) and 40 depth intervals (~3m per interval). We then use the makemaps() function to generate maps of island extents, using the projection EPSG:2100 (GGRS87/Greek Grid).

#generate interval file with a time cutoff of 20 kya (LGM), and 40 intervals (~3m per interval)
getintervals_sealvl(20,40)

#generate maps of island extents
makemaps("AegeanIslands.asc",epsg=2100)

We can visualise the effect of sea level change over the past 20,000 years by animating the maps generated by PleistoDist. Note that because we used sea level depth intervals to generate the interval file, this animation is not in chronological order.

library(magick)
library(purrr)
library(gtools)
mixedsort(list.files(path="output/raster_flat/",pattern="*.tif",full.names=T)) %>% 
  map(image_read) %>% 
  image_join() %>% 
  image_animate(fps=2) %>% 
  image_write("aegeanislands_animated.gif")

We need to generate a shapefile of sampling points for calculating island area over time. To do this, we will use the ‘Add Delimited Text Layer’ function in QGIS to read the Aegean_Angiosperms.csv file, rename the ‘Islands’ colume to ‘Name’, and filter for rows where typelgm == I (filtering out islands that were previously connected to the mainland).

We can now use the pleistoshape functions in PleistoDist to calculate the island area for each sea level interval as well as the mean island area over time. To see if elevation plays a role in affecting species richness, we can also use PleistoDist to calculate the surface area of each island.

#calculate island area over time
pleistoshape_area("Aegean_Angiosperms_islandsonly.shp",2100)
pleistoshape_surfacearea("Aegean_Angiosperms_islandsonly.shp",2100)

After calculating the mean island area and surface area, we can use simple linear models to assess the relationship between island area and species richness. We will exclude the island of Skantzoura since the sampling point fell outside the inferred island boundary for at least the first 6 intervals.

#load and filter PleistoDist outputs
islandarea <- read.csv("output/island_area.csv") %>%
  dplyr::select(Island,mean) %>%
  dplyr::filter(Island != "Skantzoura")
islandsurfacearea <- read.csv("output/island_surfacearea.csv") %>%
  dplyr::select(Island,mean) %>%
  dplyr::filter(Island != "Skantzoura")

#load species richness dataset from Hammoud et al (2021) and rename island name column for consistency
angiosperms <- read.csv("Aegean_angiosperms.csv",sep = ";") %>%
  dplyr::select(island,Sn,Nn,En)
colnames(angiosperms)[1] <- "Island"

#join PleistoDist outputs with species richness data from Hammoud et al (2021) dataset
islandarea <- left_join(islandarea,angiosperms,by="Island")
islandsurfacearea <- left_join(islandsurfacearea,angiosperms,by="Island")

#run linear models
area_allsp <- lm(log(Sn) ~ log(mean),data=islandarea)
surfacearea_allsp <- lm(log(Sn) ~ log(mean),data=islandsurfacearea)
area_native <- lm(log(Nn) ~ log(mean),data=islandarea)
surfacearea_native <- lm(log(Nn) ~ log(mean),data=islandsurfacearea)
area_endemic <- lm(log(En) ~ log(mean),data=islandarea)
surfacearea_endemic <- lm(log(En) ~ log(mean),data=islandsurfacearea)
Chorotype 2D Area Adjusted \(R^2\) 2D Area p-value 2D Area AIC Surface Area Adjusted \(R^2\) Surface Area p-value Surface Area AIC
Total Species Richness 0.3917022 7.0415557^{-7} 68.9961019 0.392509 6.8147118^{-7} 68.9297401
Native Non-endemic Species Richness 0.3873976 8.3800519^{-7} 70.1025982 0.3882024 8.1124869^{-7} 70.0368659
Endemic Species Richness 0.2662759 7.4543655^{-5} 99.7121938 0.2665386 7.3875836^{-5} 99.6942878

Linear models show that species-area relationships are significantly positively correlated for all angiosperm chorotypes across the true islands of the Aegean Sea, which is consistent with the findings of Hammoud et al (2021). However, our extended analysis also shows that linear models incorporating surface area instead of 2D area show a marginally better fit to the species richness data for all angiosperm chorotypes, based on the AIC values for each model.

areaplot <- ggplot(islandarea)+
  geom_point(aes(x=log(mean),y=log(Nn),col="Native Species"))+
  geom_smooth(aes(x=log(mean),y=log(Nn),col="Native Species"),method="lm",se=FALSE)+
  geom_point(aes(x=log(mean),y=log(En),col="Endemic Species"))+
  geom_smooth(aes(x=log(mean),y=log(En),col="Endemic Species"),method="lm",se=FALSE)+
  labs(x="log(2D Island Area)", y = "log(Species Richness)", title = "Species-area curve for Aegean Island angiosperms")

surfaceareaplot <- ggplot(islandsurfacearea)+
  geom_point(aes(x=log(mean),y=log(Nn),col="Native Species"))+
  geom_smooth(aes(x=log(mean),y=log(Nn),col="Native Species"),method="lm",se=FALSE)+
  geom_point(aes(x=log(mean),y=log(En),col="Endemic Species"))+
  geom_smooth(aes(x=log(mean),y=log(En),col="Endemic Species"),method="lm",se=FALSE)+
  labs(x="log(Island Surface Area)", y = "log(Species Richness)", title = "Species-surface area curve for Aegean Island angiosperms")

ggarrange(areaplot,surfaceareaplot,ncol=2,nrow=1)

References